//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

//#define PRINTOUT

#include "common.cpp"

#ifndef PRINTOUT
static void trace (Bitmap& bmap, Bitmap& uniq,
  coord_t& start, coord_t& begin, coord_t& end,
  size_t& bonds, size_t& ones, size_t& twos, size_t& threes,
  size_t& segs, size_t& sides, size_t& ends, size_t& lines, size_t& corners)
{
  bonds = ones = twos = threes = 0;
  segs = sides = ends = 0;
  lines = corners = 0;

  size_t mbits = 0;

  coord_t cur = begin;
  while ( 1 )
    {
      int rx = cur.pos.x - start.pos.x;
      int ry = cur.pos.y - start.pos.y;
      bool bounce = bmap.get (rx, ry);
      bool washere = uniq.get (rx, ry);

      mbits = (mbits << 1) | (bounce ? 1 : 0);
      if ( (mbits & 0x1f) == 0x11 )
          ends++;

      if ( (mbits & 0x07) == 0x05 )
          lines++;

      if ( (mbits & 0x03) == 0x03 )
          corners++;

      segs++;

      if ( bounce )
        {
          if ( !washere )
              bonds++;
          else
              sides++;
          cur.hamster();
        }
      else
        {
          if ( !washere )
            {
              int cnt = 0;
              if ( Field::upw (rx, ry) )
                {
                  cnt += bmap.get (rx, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry - 1) ? 1 : 0;
                }
              else
                {
                  cnt += bmap.get (rx, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry + 1) ? 1 : 0;
                }
              if ( cnt == 1 )
                  ones++;
              else if ( cnt == 2 )
                  twos++;
              else if ( cnt == 3 )
                  threes++;
            }
          cur.punch();
          cur.hamster();
        }
      if ( !washere )
          uniq.set (rx, ry);

      if ( cur.pos.x == end.pos.x && cur.pos.y == end.pos.y )
          break;
    }
}
#else

static inline double rotx(const vect2d& pos)
{
  return pos.x-pos.y;
}

static inline double roty(const vect2d& pos)
{
  return pos.x+pos.y;
}

static void print (const std::string& base, Bitmap& bmap, Bitmap& uniq, coord_t& start, coord_t& begin, coord_t& end)
{
  FILE *bondsfp = fopen ((base + "bonds.csv").c_str(), "w");
  FILE *endsfp = fopen ((base + "ends.csv").c_str(), "w");
  FILE *sidesfp = fopen ((base + "sides.csv").c_str(), "w");
  FILE *tracefp = fopen ((base + "trace.csv").c_str(), "w");
  FILE *linesfp = fopen ((base + "lines.csv").c_str(), "w");
  FILE *cornersfp = fopen ((base + "corners.csv").c_str(), "w");
  FILE *onesfp = fopen ((base + "ones.csv").c_str(), "w");
  FILE *twosfp = fopen ((base + "twos.csv").c_str(), "w");
  FILE *threesfp = fopen ((base + "threes.csv").c_str(), "w");
  size_t mbits = 0;
  coord_t cur = begin;
  int segstatus = 0;
  vect2d sp(0,0);
  while ( 1 )
    {
      int rx = cur.pos.x - start.pos.x;
      int ry = cur.pos.y - start.pos.y;
      bool bounce = bmap.get (rx, ry);
      bool washere = uniq.get (rx, ry);

      vect2d vpos (rx, ry);
      vect2d dir;
      if ( (cur.pos.x ^ cur.pos.y) & 1 )
          dir = vect2d (0.5, 0.5);
      else
          dir = vect2d (0.5, -0.5);

      vect2d hdir, cdir;
      if ( (cur.pos.x ^ cur.pos.y) & 1 )
        {
          if ( cur.low )
            {
              hdir = vect2d(1, 1);
              cdir = vect2d(1, -1);
            }
          else
            {
              hdir = vect2d(-1, -1);
              cdir = vect2d(-1, 1);
            }
        }
      else
        {
          if ( cur.low )
            {
              hdir = vect2d(1, -1);
              cdir = vect2d(-1, -1);
            }
          else
            {
              hdir = vect2d(-1, 1);
              cdir = vect2d(1, 1);
            }
        }

      // hull segment trace
      {
      #define HSEP 0.1
      #define TSEP 0.1
      vect2d hmid = vpos - HSEP * cdir;
      vect2d cmid = vpos + HSEP * cdir;
      vect2d pp, np;

      if ( bounce )
        {
          pp = hmid - TSEP * hdir;
          np = hmid + TSEP * hdir;
        }
      else
        {
          pp = hmid - TSEP * hdir;
          np = cmid - TSEP * hdir;
        }

      fprintf (tracefp, "%g %g\n", rotx(pp), roty(pp));
      fprintf (tracefp, "%g %g\n", rotx(np), roty(np));
      }

      mbits = (mbits << 1) | (bounce ? 1 : 0);
      if ( (mbits & 0x1f) == 0x11 ) // ends
        {
          double sign;
          if ( (cur.pos.x ^ cur.pos.y) & 1 )
              sign = (bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1) || bmap.get (rx + 1, ry)) ? -1 : 1;
          else
              sign = (bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx + 1, ry)) ? -1 : 1;
          fprintf (endsfp, "%g %g\n\n", rotx(vpos + sign * dir), roty(vpos + sign * dir));
        }

      if ( (mbits & 0x07) == 0x05 ) // lines
        {
          vect2d pp = vpos - hdir;
          vect2d np = vpos;
          fprintf (linesfp, "%g %g\n", rotx(pp), roty(pp));
          fprintf (linesfp, "%g %g\n\n", rotx(np), roty(np));
        }

      if ( (mbits & 0x03) == 0x03 ) // corners
        {
          vect2d pp = vpos - 0.5 * hdir;
          fprintf (cornersfp, "%g %g\n\n", rotx(pp), roty(pp));
        }

      if ( bounce )
        {
          // segment trace printout
          {
            vect2d pos1 = vpos + dir;
            vect2d pos2 = vpos - dir;
            if ( segstatus == 0 )
              {
                if ( cur.low )
                    sp = pos2;
                else
                    sp = pos1;
                fprintf (bondsfp, "%.2f %.2f\n", rotx(sp), roty(sp));
                segstatus++;
              }

            if ( (sp - pos1).length() > 0.01 )
                sp = pos1;
            else
                sp = pos2;

            fprintf (bondsfp, "%.2f %.2f\n", rotx(sp), roty(sp));
          }

          if ( !washere )
            {
              //bonds++;
            }
          else // sides
            {
              fprintf (sidesfp, "%g %g\n\n", rotx(vpos), roty(vpos));
            }
          cur.hamster();
        }
      else
        {
          if ( !washere )
            {
              int cnt = 0;
              if ( Field::upw (rx, ry) )
                {
                  cnt += bmap.get (rx, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry - 1) ? 1 : 0;
                }
              else
                {
                  cnt += bmap.get (rx, ry + 1) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx + 1, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx, ry - 1) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry) ? 1 : 0;
                  cnt += bmap.get (rx - 1, ry + 1) ? 1 : 0;
                }

              if ( cnt == 1 )
                  fprintf (onesfp, "%g %g\n\n", rotx(vpos), roty(vpos));
              else if ( cnt == 2 )
                  fprintf (twosfp, "%g %g\n\n", rotx(vpos), roty(vpos));
              else if ( cnt == 3 )
                  fprintf (threesfp, "%g %g\n\n", rotx(vpos), roty(vpos));
            }
          cur.punch();
          cur.hamster();
        }
      if ( !washere )
          uniq.set (rx, ry);

      if ( cur.pos.x == end.pos.x && cur.pos.y == end.pos.y )
          break;
    }

  fclose (bondsfp);
  fclose (endsfp);
  fclose (sidesfp);
  fclose (tracefp);
  fclose (linesfp);
  fclose (cornersfp);
  fclose (onesfp);
  fclose (twosfp);
  fclose (threesfp);
}
#endif

static void usage (const char *argv0)
{
  fprintf (stderr, "%s: %s [-e exponent] [-n] [-p factor] [-r seed] [-R rnd.bin] [-f] [-o perc.csv] [-s samplecount] [-i start end] [-g generator] [-b] [-l growthpow startlayer endlayer]\n", argv0, argv0);
  fprintf (stderr, "    -e exponent - set the correlation exponent\n");
  fprintf (stderr, "    -n - non-progressive calculation\n");
  fprintf (stderr, "    -p factor - specify the progression factor\n");
  fprintf (stderr, "    -r seed - set the random seed\n");
  fprintf (stderr, "    -R rnd.bin - random binary bitstream data\n");
  fprintf (stderr, "    -o perc.csv - output data to file\n");
  fprintf (stderr, "    -s samplecount - how many samples to do at first level\n");
  fprintf (stderr, "    -i start end - calculate with A from start to end\n");
  fprintf (stderr, "    -g generator - use the given random generator, put 'list' to see available\n");
  fprintf (stderr, "    -b - add baseline\n");
  fprintf (stderr, "    -m - use only 0 potential for sealevel and no perimeter edge test\n");
  fprintf (stderr, "    -l growthpow startlayer endlayer - specify the layer growth base and layer range\n");
  fprintf (stderr, "    -h - trace only half\n");
  fprintf (stderr, "    -H - trace only the hull\n");
  fprintf (stderr, "    -f - generate frequency formated data instead of the aggregate\n");
  exit(-1);
}

int main (int argc, char *argv[])
{
  big_t SAMPLES = 16384;
  bool progressive = true;
  const char *outfn = "perc.csv";
  int start = 32;
  int end = 2048;
  int rnd_seed = 0;
  EXPONENT = 0.75;
  bool add_baseline = false;
  double pf = 2.0;
  const char *rndpath = 0;
  bool middle = false;
  bool dohalf = false;
  bool onlyhull = false;
  bool freq = false;

  const gsl_rng_type * rng_type = gsl_rng_default;

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' || argv[i][1] == 0 || argv[i][2] != 0 )
          usage (argv[0]);

      switch (argv[i][1])
        {
          case 'i': {
            if ( i + 2 >= argc )
                usage(argv[0]);
            if ( sscanf (argv[++i], "%d", &start) != 1 ||
                sscanf (argv[++i], "%d", &end) != 1 ||
                start < 1 || end < start )
                usage (argv[0]);
            break;
          }
          case 'n': {
            progressive = false;
            break;
          }
          case 'b': {
            add_baseline = true;
            break;
          }
          case 'h': {
            dohalf = true;
            break;
          }
          case 'H': {
            onlyhull = true;
            break;
          }
          case 'f': {
            freq = true;
            break;
          }
          case 'm': {
            middle = true;
            break;
          }
          case 'o': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            outfn = argv[i];
            break;
          }
          case 'r': {
            i++;
            if ( i == argc || sscanf (argv[i], "%d", &rnd_seed) != 1 )
                usage (argv[0]);
            break;
          }
          case 'R': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rndpath = argv[i];
            break;
          }
          case 'p': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &pf) != 1 )
                usage (argv[0]);
            break;
          }
          case 'e': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &EXPONENT) != 1 )
                usage (argv[0]);
            break;
          }
          case 'l': {
            if ( i++ == argc || sscanf (argv[i], "%lf", &BASEPOW) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYERSTART) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYEREND) != 1 )
                usage (argv[0]);
            break;
          }
          case 's': {
            i++;
            if ( i == argc || sscanf (argv[i], "%llu", &SAMPLES) != 1 )
                usage (argv[0]);
            break;
          }
          case 'g': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rng_type = 0;
            for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
              {
                if ( strcmp ((*t)->name, argv[i]) == 0 )
                  {
                    rng_type = *t;
                    break;
                  }
              }
            if ( !rng_type )
              {
                printf ("Random generators available:\n");
                for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
                    printf ("  %s\n", (*t)->name);
                exit(0);
              }
            break;
          }
          default: {
            usage (argv[0]);
            break;
          }
        }
    }

  if ( end > FIELD_SIZE )
    {
      fprintf (stderr, "Error: interval too large for the compile time field size!\n");
      abort();
    }

  rndf.seed (rng_type, rnd_seed, rndpath);

  printf ("Using random generator: %s\n", rng_type->name);
  printf ("Random seed at %d\n", rnd_seed);
#ifdef UNCORRELATED
  printf ("Uncorrelated, ignoring the correlation exponent\n");
#else
  printf ("Correlation exponent is %f\n", EXPONENT);
  printf ("Cumulation %d\n", CUMULATION);
  printf ("Layer size growth power %g\n", (double)BASEPOW);
  printf ("Adding baseline: %s\n", add_baseline ? "YES" : "NO");
#endif
  printf ("Field size %d\n", FIELD_SIZE);


  //calc_cell_dist (outfn);
  //calc_dist (outfn);
  //exit(0);

  Field *f = new Field(add_baseline);

  //exit(0);
  //f->printsurf("surf.dat", 4096, 2048, false);
  //f->testdim();

  printf ("Running, %s from %llu samples\n", progressive ? "progressive" : "non-progressive", SAMPLES);

#ifndef PRINTOUT
  FILE *fpcsv = fopen (outfn, "w");
  if ( !fpcsv )
      exit(-1);

  printf ("Writing data into '%s'..\n", outfn);

  if ( freq )
      fprintf (fpcsv, "A,PROPERTY,VALUE,COUNT\n");
  else
      fprintf (fpcsv, "A(BONDS,FULLS,LINES),M,SAMPLES,SUM(V1),SUM(V2),...,SUM(VM),"
          "SUM(V1*V1),SUM(V1*V2),...,SUM(V1*VM),SUM(V2*V2),SUM(V2*V3),...,SUM(VM*VM)\n");
  fflush (fpcsv);
#endif

  aggr_freq aggrf;

  double SC = SAMPLES;

  for ( double scale = start; scale <= end; scale *= pf, SC /= (progressive ? pf : 1.0) )
    {
      SAMPLES = SC;
      int A = scale;

#ifndef PRINTOUT
      printf ("Calculating A=%d, %llu samples..\n", A, SAMPLES);
#endif

      int first = 0;
      int last = first + A + 2;
      int mid = (first + last) / 2;
      int offset = FIELD_SIZE / 2 - mid;

      aggregator<18> aggr_all;
      aggregator<9> aggr_hull;
      size_t samples = 0;
      while ( SAMPLES > samples )
        {
          f->reset();
          int x = mid;
          int y = mid;

          coord_t start;
          start.pos.x = offset + x;
          start.pos.y = offset + y;
          start.low = rndf.rnd() >= 0;

          if ( !middle )
            {
#ifndef UNCORRELATED
              f->set_sealevel (f->value (offset + x, offset + y));
#else
              if ( !f->valueb (offset + x, offset + y) )
                  continue;
#endif
            }
          else
            {
              if ( !f->valueb (offset + x, offset + y) )
                  continue;
            }

          bool left, right, up, down;
          left = right = up = down = false;

          // setup
          int limit = last - first;
          static Bitmap bmap(32); // hull bond map
          static Bitmap bmap2(32); // perimeter bond map
          static Bitmap bmap3(32); // hull trace unique bond map
          static Bitmap bmap4(32); // perimeter trace unique bond map
          bmap.reset (limit + 1);
          bmap2.reset (limit + 1);
          bmap3.reset (limit + 1);
          bmap4.reset (limit + 1);
          int hlimit = limit / 2;
          coord_t save;

          // first wind forward
          coord_t cur = start;
          bool looped = false;
          while ( 1 )
            {
              if ( start.pos.x - cur.pos.x == hlimit - 1 || cur.pos.x - start.pos.x == hlimit ||
                  start.pos.y - cur.pos.y == hlimit - 1 || cur.pos.y - start.pos.y == hlimit )
                  break;
              //if ( abs (start.pos.x - cur.pos.x) == hlimit || abs (start.pos.y - cur.pos.y) == hlimit )

              save = cur;

              if ( f->valueb (cur.pos.x, cur.pos.y) )
                {
                  bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
                  cur.hamster();
                }
              else
                {
                  cur.punch();
                  cur.hamster();
                }
              if ( cur == start )
                {
                  looped = true;
                  break;
                }
            }

          if ( looped )
            {
              //printf ("Looped!\n");
              continue;
            }

          coord_t end = save;
          if ( !bmap.get (end.pos.x - start.pos.x, end.pos.y - start.pos.y) )
              end.punch();
          end.hamster();

          if ( dohalf )
            {
              end = start;
              end.hamster();
              bmap.reset(limit + 1);
            }

          // trace backwards too
          cur = start;
          while ( 1 )
            {
              if ( start.pos.x - cur.pos.x == hlimit - 1 || cur.pos.x - start.pos.x == hlimit ||
                  start.pos.y - cur.pos.y == hlimit - 1 || cur.pos.y - start.pos.y == hlimit )
                  break;
              //if ( abs (start.pos.x - cur.pos.x) == hlimit || abs (start.pos.y - cur.pos.y) == hlimit )

              save = cur;

              if ( f->valueb (cur.pos.x, cur.pos.y) )
                {
                  bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
                  cur.rhamster();
                }
              else
                {
                  cur.punch();
                  cur.rhamster();
                }
            }

          coord_t begin = save;
          if ( !bmap.get (begin.pos.x - start.pos.x, begin.pos.y - start.pos.y) )
              begin.punch();

          if ( !onlyhull )
            {
              bool pfail = false;
              // tracing the perimeter
              cur = begin;
              while ( 1 )
                {
                  int rx = cur.pos.x - start.pos.x;
                  int ry = cur.pos.y - start.pos.y;
                  bool bounce = bmap.get (rx, ry);

                  //printf ("Looping: %d %d %d\n", rx, ry, bounce);
                  if ( !bounce )
                    {
                      if ( (cur.pos.x ^ cur.pos.y) & 1 )
                        {
                          if ( cur.low )
                            {
                              if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1)) &&
                                  (bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx - 1, ry)))
                                  bounce = true;
                            }
                          else
                            {
                              if ((bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1) || bmap.get (rx + 1, ry)) &&
                                  (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1)))
                                  bounce = true;
                            }
                        }
                      else
                        {
                          if ( cur.low )
                            {
                              if ((bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx + 1, ry)) &&
                                  (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry + 1) || bmap.get (rx, ry + 1)))
                                  bounce = true;
                            }
                          else
                            {
                              if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1)) &&
                                  (bmap.get (rx, ry + 1) || bmap.get (rx - 1, ry + 1) || bmap.get (rx - 1, ry)))
                                  bounce = true;
                            }
                        }
                    }

                  if ( bounce )
                    {
                      bmap2.set (rx, ry);
                      cur.hamster();
                    }
                  else
                    {
                      cur.punch();
                      cur.hamster();
                    }

                  if ( cur.pos.x == end.pos.x && cur.pos.y == end.pos.y )
                      break;

                  if (cur == begin)
                    {
                      pfail = true;
                      break;
                    }
                }

              if ( pfail )
                {
                  //printf ("Pfail continue\n");
                  continue;
                }
            }

#ifdef PRINTOUT
          print ("hull_", bmap, bmap3, start, begin, end);
          if ( !onlyhull )
              print ("perim_", bmap2, bmap4, start, begin, end);
#else
          size_t hull_bonds, hull_ones, hull_twos, hull_threes, hull_segs, hull_sides, hull_ends, hull_lines, hull_corners;
          trace (bmap, bmap3,
              start, begin, end,
              hull_bonds, hull_ones, hull_twos, hull_threes,
              hull_segs, hull_sides, hull_ends, hull_lines, hull_corners);

          size_t perim_bonds, perim_ones, perim_twos, perim_threes, perim_segs, perim_sides, perim_ends, perim_lines, perim_corners;
          if ( !onlyhull )
          trace (bmap2, bmap4,
              start, begin, end,
              perim_bonds, perim_ones, perim_twos, perim_threes,
              perim_segs, perim_sides, perim_ends, perim_lines, perim_corners);

          if ( freq )
            {
              if ( onlyhull )
                {
                  aggrf.add (
                    A,
                    hull_bonds, hull_ones, hull_twos, hull_threes, hull_segs, hull_sides, hull_ends, hull_lines, hull_corners
                  );
                }
              else
                {
                  aggrf.add (
                    A,
                    hull_bonds, hull_ones, hull_twos, hull_threes, hull_segs, hull_sides, hull_ends, hull_lines, hull_corners,
                    perim_bonds, perim_ones, perim_twos, perim_threes, perim_segs, perim_sides, perim_ends, perim_lines, perim_corners
                  );
                }
            }
          else
            {
              if ( onlyhull )
                {
                  aggr_hull.add (
                    hull_bonds, hull_ones, hull_twos, hull_threes, hull_segs, hull_sides, hull_ends, hull_lines, hull_corners
                  );
                }
              else
                {
                  aggr_all.add (
                    hull_bonds, hull_ones, hull_twos, hull_threes, hull_segs, hull_sides, hull_ends, hull_lines, hull_corners,
                    perim_bonds, perim_ones, perim_twos, perim_threes, perim_segs, perim_sides, perim_ends, perim_lines, perim_corners
                  );
                }
            }
#endif

          //printf ("Found: %zu %d %d %d %d\n", cluster_size, left, right, up, down);
#ifdef PRINTOUT
          break;
#endif
          samples++;
        }
#ifndef PRINTOUT
      assert ( samples == SAMPLES);
      if ( !freq )
        {
          if ( onlyhull )
              fprintf (fpcsv, "%d,%s\n", A, aggr_hull.sprint());
          else
              fprintf (fpcsv, "%d,%s\n", A, aggr_all.sprint());
        }
      fflush (fpcsv);
#else
      break;
#endif
    }
  delete f;

#ifndef PRINTOUT
  if ( freq )
      aggrf.print(fpcsv);

  fclose (fpcsv);
#endif
  return 0;
}

